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In the first part of this paper, a review is given on the mechanism for the disappearance 
of an intrinsic stacking fault in a hard-sphere (HS) crystal under gravity, which we recently 
discovered by Monte Carlo (MC) simulations [A. Mori et at, J. Chem. Phys. 124 (2006), 
17450; Mol. Phys. 105 (2007), 1377]. We have observed, in the case of fee (001) stacking, 
that the intrinsic stacking fault running along an oblique direction shrunk through the gliding 
of a Shockley partial dislocation at the lower end of the stacking fault. In order to address 
the shortcomings and approximations of previous simulations, such as the use of periodic 
boundary condition (PBC) and the fact that the fee (001) stacking had been realized by 
the stress from the small PBC box, we present an elastic strain energy calculation for an 
infinite system and a MC simulation result for HSs in a pyramidal pit under gravity. In 
particular, the geometry of the latter has already been tested experimentally [S. Matsuo 
et al., Appl. Phys. Lett. 82 (2003), 4283]. The advantage of using a pyramidal pit as a 
template as well as the feasibility of the mechanism we describe is demonstrated. 

§1. Introduction 

In 1957, the existence of the crystalline phase in a hard-sphere (HS) system 
was found by Monte Carlo (MC) and molecular dynamics (MD) simulations the 
crystalline phase transition in the HS system is sometimes referred to as the Alder 
transition. It was surprising that a system comprised merely of hard-core repulsion 
exhibited a phase transition. An intuitive understanding of the Alder transition is 
given by decomposing of the entropy into a contribution due to the configurational 
variation of the HSs' centers and a contribution due to the vibrational degree of 
freedom around the HSs' equilibrium positions. While the configurational entropy 
dominates in a disordered fluid phase, the vibrational entropy dominates in the 
crystalline phase. In the HS system the phase behavior is governed by the density; 
the system is in the fluid phase at a density lower than (j) f and in the crystalline phase 
of a face-centered cubic (fee) structure at a density higher than <p s . Here, the particle 
density is expressed in terms of the volume fraction of the HSs <fi = na 3 N/QV, where 
a is the HS diameter, N is the number of particles, and V is the total system volume. 
A crystal of <p s can coexist with a fluid of <j)f when the volume fraction of the total 
system lies between cftf and <f> a . The Hoover and Ree^ 1 determined </>/ = 0.494 and 
<j) s = 0.545 by a MC method in 1968. Those values have been revised in the last 
decade to <pf = 0.491 and <p s = 0.542 by a MD simulation study of the direct two- 
phase coexistence (i.e., the crystal/fluid interface).® We note that even the first MD 
simulation of the HS crystal/fluid interface was successfully performed within the 
last decade P 
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The present situation is a little different from that in the early years. In the 
1960s and -70's, the existence of colloidal crystals drew much interest as an experi- 
mental realization of the Alder transition. The effective HS picture was proposed for 
charged colloids, which interacts through a screened Coulomb potential where the 
interparticle interaction is thus well described by a repulsive Yukawa form.^ Today, 
the HS system is not just an effective model of the colloids; poly(methylmethacrylate) 
(PMMA) particles with stabilizing polymers grafted on the particle's surface® were 
developed as HS suspensions. The PMMA particles were dispersed in a compounded 
hydrocarbon medium so that the HS nature with regards to the crystal-melt phase 
transition was exhibitedP'^'^ffl'EQ 

It should also be noted that colloidal crystals can possibly be used as materials 
for photonic crystals.^ To realize a photonic band, the defects in the colloidal 
crystal should be reduced. To this end, many techniques have been developed. One 
of them is colloidal epitaxy, in which a patterned substrate is used as a template to 
fix the stacking direction to be [001] in the sedimentation of the colloidal particlesPSJ 
Various patterns for the template have been examined, but in this work only a single 
pyramidal pitP^ is considered. 

In order to improve the quality of the colloidal crystal, the effect of gravity on 
the crystallinity of the colloidal crystal must be taken into account. Zhu et a/P*^ 
performed a colloidal crystallization on the Space Shuttle and concluded that under 
microgravity, a random hexagonal close pack (rhcp) was found. On the oher hand, 
under normal gravity, the colloidal crystal formed by sedimentation is a rhcp/fcc 
mixture although the colloids sometimes freeze into a glassy state. Here, rhcp is 
the random staking of hexagonal planes (fee {111} or hep (0001)); viewed along 
(111), the fee is of ABCABC- • • type, and the hep is of ABAB- • • type while the 
rhcp corresponds to a random sequence. In other words, the stacking disorder can 
be reduced by gravity. We have, however, not found a final answer to the stacking 
sequence of colloidal crystals under gravity. Although the gravitational constant 
g* = mga/kBT was different, Kegel and DhonlP^ observed a faulted twinned fee 
under gravity. Here, m is the particle's mass, g is the acceleration due to gravity, a 
is the HS diameter, and fcgT is the temperature multiplied by Boltzmann's constant. 
We note that the temperature T should be defined from the thermal motion of the 
dispersing particles. 

This paper focuses on the mechanism of the disappearance of the stacking disor- 
der in a HS crystal under gravity. We have already demonstrated the disappearance 
of stacking disorder in a HS crystal under gravity by MC simulations P® We have 
examined closely where the shrinking of an intrinsic stacking fault was observed.^ 
In §[21 we review those simulation results, point out their shortcomings, and reiterate 
the mechanism of the defect disappearance that had been found. The objective of 
this paper is to provide new results to complement these shortcomings. The sys- 
tems^ 'D^ were highly stressed because of their smallness. We preset some results 
of the elastic energy calculation for an infinite system in § [3j The periodic boundary 
condition (PBC) applied in the previous simulations^'^ was one of the artifacts 
that needed to be addressed. In § 0] we present results of the MC simulation for 
realizable boundary conditions by simulating HSs in a pyramidal pitP^ 
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§2. Monte carlo simulation under the PBC 

In previous simulations we pre- 
sented results for two system sizes. One 
contained A r =1664 particles and the 
other 7V=3744. In the former, fee (001) 
stacking occurred, while in the latter, 
stacking of the hexagonal layers took 
place. We concentrated on the for- 
mei HEli Decause the defects' disappear- 
ance was observed in this case. Al- 
though the side length equaled 4ao, six 
particles lay along a side. Here, ao = 
1.57o- = (4<7 3 /Vs) 1/3 was the fee lat- 
tice constant at the crystal/fluid coexis- 
tence of Hoover and Ree's value® with 
p denoting the particle density [<fi = 
(7r/6)/9<7 3 ]. Six unit cells lay along the 
diagonal direction. In other words, the 
[010] and [100] are directed to the diago- 
nal directions P^* This indicated that the 
crystal at the bottom was at a higher 
pressure than the crystal-fluid coexis- 
tence pressure [the Hoover and Ree's 
value is (Pa 3 /k B T) coex = 11.7P and 
Davidchack and Laird's value 11.5 5^^1 in 
accordance with the mechanical balance 
equation 

dP/dz 

at the sedimentation equilibrium, where p(z) is the particle density on a coarse-scale 
at altitude z. Also, (|2-ip implies that the lattice spacing along z direction varies 
and that the unit cell is compressed along the vertical direction, i.e., the crystal 
is no longer a cubic system, as was previously confirmed.^ 1 The disappearance of 
the stacking fault occurred under such artificially stressed condition. Such stress, 
however, plays a central role in defect reduction. As shall be described below, the 
(001) stacking, in which the shrinking of the stacking fault occurred,^ was thereby 
induced. Moreover, some researchers made an efforts to realize the (001) growth, 
such as the colloidal epitaxy,^ with the understanding that the stacking sequence 
for this growth is unique. Concentrating on phenomena occurring under the (001) 
stacking without particularly concerning ourselves with the experimental realization 
of this stacking is one of paths to elucidation of the defect reduction mechanism. 

In previous simulations E^ED the gravitational constant, g*, was changed in a 
step- wise fashion in order to avoid trapping the system in a metastable state such as 
a polycrystalline stateP^ We set g* to a certain value for 2xl0 5 MC cycle (MCC), 
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and then changed g* by an amount of Ag* = 0.1, where one MCC is defined so as to 
include, on average, one position move per particle. Transformation of a defective 
crystal into a less-defective crystal was observed for the case of (001) stacking for 
g* around g* = 0.9.G3 While g* was kept at 0.9 for the (001) stacking case, we 
found that the shrinking of an intrinsic stacking fault was mediated by the glide of 
the Shockley partial dislocation terminating at the bottom end of the stacking fault 
running in the oblique direction.^ In Fig. [T]we see a disappearance of the stacking 
faults, which run from the lower-left to the middle-right, [see, Ref. 18) for details of 
the shrinking of the intrinsic stacking fault]. We stress that the advantage of the 
(001) stacking is not only the uniqueness of the stacking sequence but also the glide 
mechanism of the Shockley partial dislocation. 



§3. Elastic strain energy consideration 



In this section, we present the elas- 
tic energy calculation for a system in- 
cluding an intrinsic stacking fault run- 
ning along the [111]. The lower end 
of the stacking fault is terminated by a 
Shockley partial dislocation, such as the 
one shown in Ref. 18) and illustrated in 
Fig. [2j In this paper, we incorporated 
the effect of gravity as buoyancy due to 
the particle deficiency accompanied by 
the dislocation core. The particle defi- 
ciency is one-third in Fig. [2j Though 
the crystal was indeed strained due to 
gravity,^ 1 which coupled with the stress 
due to the Shockley partial dislocation 
to help promote the shrinking of the in- 
trinsic stacking fault, we can understand 
the shrinking of the intrinsic stacking fault 
between these stresses. 

The Shockley partial dislocation 
terminating an intrinsic stacking fault 
running along the [111] (as shown in 
Fig. [2]) is defined by the Burgers vector 
b 1 = (1/6)[211] = -ai/3 + a 2 /6 + a 3 /6 
with oi, ot2, and a 3 being the lattice 
vectors, which is shown by the arrow 
in Fig. [3l We can understand the par- 
tial dislocation by the decomposition of 
a perfect dislocation. In Fig. O the 
Burgers vector of a perfect dislocation, 
b = (1/2) [110], is decomposed as b = 



1] Shockley partial viewed from [110] 

mm 



Fig. 2. Illustration of an intrinsic stacking 
fault. Particles in the distorted crystal 
are indicated by dots, and the regular lat- 
tice positions are indicated by open circles. 
The dotted line indicates the stacking fault. 
For simplicity, particles outside the portion 
right of stacking fault are not displaced in 
this illustration. 

without considering the cross coupling 



[001 



010] 




; #-[110] 



[100] 



Fig. 3. The Burgers vector b 1 = (1/6) [211] 
(arrow) and the (111) plane (painted). The 
arrow connects a lattice position to an ad- 
jacent lattice position, say lattice point A 
to lattice point B. 
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b 1 + b 11 with b 11 = (1/6)[121] as was done in Ref. 151]) . 

We calculated the elastic energy due to a dislocation running along a unit vector 
£ using the following formula from isotropic, linear elastic theory® 

„,,„. ub 2 ( ~ sin 2 #\ ( aR\ 

where p,, u, and 6 are the shear modulus, the Poisson ratio, and the angle between 
the Burgers vector b and £, respectively, and b denotes \b\. Here, with a, which is 
several tenths, the radius of defined core region is ro = b/a, and R is the dimension 
over which the elastic field expands, which we can usually identify with the crystallite 
or grain radius. In the present case, where the stacking fault is running upward along 
the (111) starting at the position of the Shockley partial dislocation, we simply set R 
to the distance from the upper boundary to the Shockley partial dislocation. (Note 
that the dependence of the geometry of the boundary or the shape of the crystallite 
is ignored.) Substituting b 1 for b (i.e., |b 7 | = a/y/6 with a being the fee lattice 
constant) and £ = (l/v2)[110] (i.e., 9 = tt/6), we obtain an elastic energy of 

Hereafter, we consider a system of unit length thickness perpendicular to Fig.[2j The 
core region is defined so that linear elastic theory is still valid outside that region. 
Borrowing the empirical result for metals,® the core energy U core is proportional to 
[ib 2 . Let us write U core = j3pb 2 with j3 being a certain constant of order less than 
unity. We have 

U core = f3pa 2 /6, (3-3) 

(b = b 1 ) for the intrinsic stacking fault. Note that ()3-3[) is independent of R. 

We should now consider the stacking fault energy. This quantity, as well as the 
elastic and core energies (rigorously, the free energies), is an entropic contribution 
resulting from the variation of the vibrational mode distribution. Though the core 
energy has not been calculated, the shear modulus p was calculated by a MC sim- 
ulation 22 ' and density functional theory,® and the stacking fault energy j s f was 
calculated by a MC simulation.® The quantity, pa^/ksT for the HS crystal ranges 
between 50 and 100, depending on the particle number density, p. This range of 
p corresponds to pa 3 = 1.06 — 1.13 (a/a = 1.55 - 1.52) where the disappearance 
of the stacking disorder was observed in the MC simulations.^'^ The stacking 
fault interfacial energy per unit area, 7s/c 2 /ksl ', was (at most) of the order 10 -4 at 
pa 3 = 1.10. The difference of the orders in p, and j s f has a crucially important mean- 
ing; if the perfect dislocation is decomposed into two partial dislocations connected 
by a stacking fault, the separation between two partial dislocations would be much 
longer than the order of the crystallite radius, which is typically a few to several 
hundred lattice spacings. The total stacking fault energy is obtained by multiplying 
the length of the stacking fault, which is C-R; the proportional constant C depends 
on the geometry of the boundary. Thus, we have U s f = ClsfR- 
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The gravitational energy is solely given by U g {R) — U g (0) = mpa 2 R/3, where 
we have neglected the dependence of the particle density on the altitude. This 
approximate treatment was consistent with the treatment where we neglected the 
deformation according to (|2- 1 j) . Here U g (R) is the gravitational energy where the 
lower end of the stacking fault is located from the upper boundary of the crystallite; 
so U g (0) vanished because the stacking fault went out of the crystallite when R = 0. 
The sum of gravitational energy and the stacking fault energy was thus a linear 
function 

U. f + U g = (Clsf + mpa 2 /3)R. (3-4) 



The total energy, U e i + Uoore + U s f + U g , is comprised of a logarithmic term (|3-2|) , 
a constant (|3-3p . and a linear term (|3-4p . The coefficient of the logarithmic term is 
positive and that of the linear term is also apparently positive. This means that as 
R decreases, the energy decreases. We have shown by an elastic calculation that in 
the (001) stacking, gravity provided a driving force that promoted the glide of the 
Shockley partial dislocation upwards. 



§4. Monte carlo simulation of hard spheres in a pyramidal pit 



In the preceding two sections, we 
considered two systems with artificial 
boundary conditions. In contrast, in 
this section, we treat a realistic bound- 
ary condition. Mastsuo et alM- used 
a pyramidal pit and groove to fix the 
stacking sequence. A pyramidal pit can 
be made by the anisotropic etching of 
the Si (001) surface (Fig.©. 

Contrary to the simulations with 
PBC, the crystallinity was less sensitive 
to how g* was controlled. This may be 
due to wetting on the edges of the pyra- 
midal pit (Fig. 03(a) ). In Fig. [5j snap- 
shots of MC simulations at g* = 0.1 and 
1.5 are shown; though At = 2 and 4 x 10 5 MCC, where g* was kept constant, were 
also tested, only the results for At = 1 x 10 5 MCC are shown. Figure \Mjo) shows that 
despite the edge wetting, the particles at the bottom were not crystallized. On the 
other hand, the particles were crystallized at the bottom of the pit, particularly, in 
the fee (001) stacking. In the pyramidal pit made of the {111} faces, which can be 
experimentally generated, the (001) stacking has been confirmed, as was already ex- 
perimentally demonstrated.^ Unfortunately, looking over the snapshots, we could 
not find the stacking faults, and thus, the glide mechanism of the stacking fault 
disappearance was not observed. In our opinion, however, this was an indication of 
this method's robustness against stacking disorder; the epitaxial growth at the late 
stage started from the edges of the pit where the wetting occurred. In addition, 



le.ekv X3.88K IB. a* 



Fig. 4. A top view of a pyramidal pit made of 
four Si {111} faces (scanning microscopy 
image) . 
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the geometry of the boundary was advantageous; even if a crystallite mismatched its 
crystallographic orientation or lateral position with the formed crystal as it nucleates 
at the bottom, the crystallite could slide along the wall. In other words, the simu- 
lation results did not rule out the glide mechanism of a stacking fault but instead 
emphazised its role. 




Fig. 5. Snapshots of MC simulation of the HSs in a pyramidal pit. (a) An overview at g*=0.1, 
where edge wetting is seen, (b) a magnification of (a), and (c) a magnified snapshot at g*=1.5, 
where crystalline order was observed. 



§5. Discussions 

Here we compare two simulations. Under PBC, the configuration is equivalent 
to the dislocation array of the same Burgers vector. On the other hand, owing to 
the concept of the mirror image, the pyramidal pit system is equivalent to an array 
of opposite charges. Whereas the former is less mobile, the dislocation in the latter 
more easily disappears, which can be understood by pair annihilation of opposite 
charges. The glide mechanism is likely implicated because it was observed in the 
former configuration. 

The cross term between the gravity-induced elastic field and the stress due to 
dislocation has not been treated in this paper. From (|2-ip . the particle density 
decreases as the altitude increases. This means that the buoyancy of the dislocation 
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core is promoted as the core goes up because the higher the altitude is, the looser 
the packing is. Essentially the same argument can be valid for the strain energy to 
due the dislocation. 

§6. Concluding remarks 



We presented two simulation results and one theoretical calculation. All of our 
results support the glide mechanism for the disappearance of the stacking disorder 
in a hard-sphere crystal under gravity. We wish to emphasize that although this 
mechanism is not the final unique answer, it is likely to occur. 

In a theoretical elastic calculation, we can understand the glide mechanism with- 
out taking into account the cross term between the gravity-induced elastic field and 
the stress due to the dislocation. That is, do~ij/dxi = has been solved, where <7jj 
is the stress tensor. The rigorous treatment is to solve d<Jij/dxi + fi = with /j 
denoting the external force (gravitational force in the present case) . This calculation 
is currently underway using the results from Ref . [T9]) . 
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Errata [Prog. Theor. Phys. Suppl. 178 (2009), 33.] 

Both in Eq. (3.2) and on the line 5 on page 38 ([6j:h page), we should add a 
prefactor (1/2^2) on mpa 2 /3 and mpa 2 R/3, respectively. This factor arizes from 
the area per [110] line on (110) plane, which is a /2V2. We can also drive this by 
considering the number of particle on [110] per unit length, i.e, U g = m(\/2/a)/3 = 
(l/2\/2)m/ja 2 /3 because p = 4/a 3 for fee lattice. 



